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Abstract 


An isothermal single-phase 3D/1D model for liquid-feed direct methanol fuel cells (DMFC) is presented. Three-dimensional (3D) mass, 
momentum and species transport in the anode channels and gas diffusion layer is modeled using a commercial, finite-volume based, compu- 
tational fluid dynamics (CFD) software complemented with user supplied subroutines. The 3D model is locally coupled to a one-dimensional 
(1D) model accounting for the electrochemical reactions in both the anode and the cathode, which provides a physically sound boundary con- 
dition for the velocity and methanol concentration fields at the anode gas diffusion layer/catalyst interface. The 1D model — comprising the 
membrane-electrode assembly, cathode gas diffusion layer, and cathode channel — assumes non-Tafel kinetics to describe the complex kinetics 
of the multi-step methanol oxidation reaction at the anode, and accounts for the mixed potential associated with methanol crossover, induced 
both by diffusion and electro-osmotic drag. Polarization curves computed for various methanol feed concentrations, temperatures, and methanol 
feed velocities show good agreement with recent experimental results. The spatial distribution of methanol in the anode channels, together with 
the distributions of current density, methanol crossover and fuel utilization at the anode catalyst layer, are also presented for different opperating 


conditions. 
© 2007 Published by Elsevier B.V. 
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1. Introduction 


Fuel cells are electrochemical devices that convert the chem- 
ical energy of an energy carrier — typically hydrogen — and 
an oxidizer — typically the oxygen of the air — directly into 
electricity and heat. In contrast to the most common pro- 
ton exchange membrane fuel cells (PEMFCs), operating with 
hydrogen, liquid-feed direct methanol fuel cells (DMFCs) use 
methanol as energy carrier, which makes them good candi- 
dates as small autonomous power sources. In fact, due to the 
high energy density of methanol, up to 100 times higher than 
state-of-the-art lithium-ion batteries, DMFCs are regarded as a 
potential substitute to conventional power generating equipment 
for portable electronic devices. 

Nevertheless, DMFCs suffer from two fundamental prob- 
lems: (i) the slow kinetics of the methanol electro-oxidation 
reaction and (11) the ability of methanol to permeate through the 
polymer membrane crossing from anode to cathode (methanol 
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crossover). In addition, there are several system design issues, 
such as water/gas/heat management or flow-field design and 
optimization, that still need a better understanding. The above 
difficulties, together with additional technological problems 
concerning auxiliar devices, such as pumps, fuel storage tanks, 
power conditioning devices, etc. have motivated a large body 
of work in the field during the last decade, combining math- 
ematical and numerical modeling with detailed experimental 
research. Particularly, the progress in DMFC modeling has been 
significant. 

Several mathematical models for DMFCs can be found in 
the literature, including early one-dimensional models [1-9] 
and more recent two- and three-dimensional models [10-14]. 
Of particular relevance is the development of models which 
account for multiphase flow and transport phenomena, such 
as the seminal work of Wang and Wang [15] or other similar 
models [16,17]. However, single-phase models such as the one 
presented here have also contributed to the understanding of 
the complex phenomena involved in DMFCs (see, e.g., Refs. 
(3-5,8,12—14,18—20]). An extensive review of this work can 
be found elsewhere [21,22] and will not be repeated here for 
brevity. 
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Nomenclature 


effective catalyst surface area per unit volume 
(m7!) 

surface area (m”) 

molar concentration (mol m~?) 

mass diffusivity of species i (m? s7!) 
electromotive force (V) 

Faraday’s constant (C mol!) 

convective mass transfer coefficient (ms~!) 
current density (A m~?) 

current (A) 

permeability (m7) 

channel length (m) 

electro-osmotic drag coefficient of species i 
molar flux (mol m~? s7!) 

pressure (Pa) 

universal gas constant (J mol`! K7!) 
source term appearing in the momentum equation 
temperature (K) 

(superficial) velocity vector (m s7 1) 

actual cell voltage (V) 

molecular weight (kg mol~!) 

coordinate along the channel (m) 
coordinate across the channel height (m) 
coordinate across the channel width (m) 


Greek letters 


overall mass transfer coefficient (m s7!) 


a 
ô thickness (m) 

€ porosity 

n overpotential (V) 

K experimental constant 

À experimental constant (mol m7?) 
u dynamic viscosity (kg m7! s7!) 
p density (kg m~?) 

o electric conductivity (S m7!) 
Subscripts 

a anode 

ac anode chanel 

acl anode catalyst layer 

agdl anode gas diffusion layer 

amb ambient 

avg average value 

c cathode 

cc cathode chanel 

ccl cathode catalyst layer 

cgdl cathode gas diffusion layer 
CO2 carbon dioxide 

cross crossover 

d drag 

e electronic 

i species i 

in channel inlet 

m methanol 


mem membrane 

O2 oxygen 

p parasitic 

ref reference value 
Ww water 
Superscripts 

0 standard conditions 
eff effective value 
m methanol 

O2 oxygen 

W water 


Typically, the anode of a liquid-feed DMFC is supplied with 
a diluted methanol aqueous solutions, while the cathode is feed 
with an oxidizer stream (air or pure oxygen) which can be either 
forced by an external blower or driven by natural convection. 
Due to the combined effect of both convection and diffusion the 
methanol and the oxygen reach the anode and cathode catalyst 
layers, respectively, where they undergo the overall electrochem- 
ical reactions: 


CH30H + H20 —> CO) + 6H% + 6e7 (1) 
cathode: 302 + 6H* + 6e7 —> 3H20 (2) 


anode: 


with standard reduction potentials Æ? = 0.02V and E? = 
1.23 V vs. saturated hydrogen electrode (SHE), respectively. 
Due to methanol crossover, reaction (1) accounts indeed for 
methanol oxidation at both the anode and the cathode, while 
reaction (2) accounts for oxygen reduction at the cathode. The 
protons generated at the anode by reaction (1) diffuse across the 
polymer membrane, while the electrons pass as current through 
the external circuit to reach the cathode, where oxygen is reduced 
with the protons and the electrons to form water according to 
reaction (2). Globally, the two electrochemical reactions are 
combined to give the overall cell reaction: 


CH30H + 302 — CO2 + 2H20 (3) 
with standard cell potential Æ? ı = 1.21 V at 298K. 


Among the performance controlling components of the direct 
methanol fuel cell, the anode electrode is known to be one of 
the most influential. In particular, the low activity of the electro- 
oxidation reaction of methanol is affected by the poisoning of 
the anode catalyst by stable adsorbed intermediates of methanol 
oxidation, which leads to large anodic overpotentials. Thus, in 
this paper we shall focus our attention on the anode side of the 
cell. 

Most published papers on DMFC modeling consider 
Butler—-Volmer or Tafel kinetics for the anodic reaction. How- 
ever, the methanol electro-oxidation reaction (1) is known to be 
a multi-step reaction that takes place as several (possibly simul- 
taneous) elementary steps at the molecular level, e.g., see Refs. 
[11,23] and references therein. 

In contrast to most early DMFC models, Meyers and New- 
man [19] proposed a mechanistic model to describe complete 
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non-Tafel methanol oxidation kinetics on Pt-Ru catalysts based 
on the sequence of elementary steps proposed by Gasteiger et 
al. [24]. An interesting feature of their model is that it appro- 
priately describes the transition from first-order kinetics for low 
methanol concentrations (high overpotentials) to zeroth-order 
kinetics for high methanol concentrations (low overpotentials), 
a transition first suggested by Ren et al. [25] and recently con- 
firmed experimentally by Vidaković et al. [23]. 

Similar kinetic models have also been used by Garcia et 
al. [8] and Kulikovsky [26] as improved models for the elec- 
trochemical oxidation of methanol. In all cases, the proposed 
kinetic expressions were used in the context of fully 1D, single- 
phase, DMFC models. Although multi-phase flow effects play 
an important role in the operation of DMFCs [15], in this paper 
the flow in the anode channels and gas diffusion layer is also 
described using the single-phase approximation, which leads to 
a simplified analytical description of the problem. Therefore, the 
aim of this paper is not to develop a mathematical model that 
fully describes all the phenomena that occur in a DMFC, but 
to illustrate the use a hybrid 3D/1D model that accounts for the 
electrochemical kinetics of Meyers and Newman [19]. 

Three-dimensional models constitute an excellent tool for 
exploring the spatial distribution of reactants and current density 
in the fuel cell. Thus, besides the overall cell polarization curves, 
we shall also study the spatial distribution of methanol in the 
anode channels, as well as the local distributions of current den- 
sity, methanol crossover and fuel utilization at the anode catalyst 
layer, of interest for the optimized design of future-generation 
DMFCs. 

The structure of the paper is as follows. The mathematical 
model is presented in Section 2; this comprises model assump- 
tions (Section 2.1), description of the physical domain (Section 
2.2), formulation of the 3D model (Section 2.3) and formu- 
lation of the 1D model (Section 2.4). The solution procedure 
is explained in Section 3, and the numerical results are pre- 
sented in Section 4, including the validation of the model against 
experiments. Finally, the concluding remarks are presented in 
Section 5. 


2. Mathematical model 
2.1. Model assumptions 


In the development of the mathematical model, a number 
of simplifying assumptions have been made. Some of them are 
central to our model, namely: 


e Carbon dioxide produced at the anode is considered to be 
dilute enough to remain dissolved in the liquid phase, i.e. 
single-phase flow is considered in the anode. 

e A large stoichiometric excess of air (or oxygen) is assumed 
in the cathode. 


Strictly speaking, the first assumption restricts the validity of 
the model to the low current density regime, when the production 
of carbon dioxide at the anode is small. However, the single- 
phase model is expected to give correct parametric trends also 


in the high current density regime. A detailed study of two- 
phase flow and transport in DMFCs results in a much more 
complex modeling approach [15], which is not the aim of this 
paper. On the other hand, the second assumption, discussed in 
Section 2.4.1 below, allows us to impose a constant value for the 
oxygen ambient concentration along the cathode channel/flow 
distributor, further simplifying the analysis. 

The rest of the assumptions, widely used in modeling studies 
of DMFCs, are the following: 


The flow is laminar and steady. 

e The temperature is constant throughout the cell. 

e The reactant concentrations are constant across the anode and 
cathode catalyst layers. 

e The concentration of methanol is sufficiently small in the 
anode for the liquid phase to be a diluted methanol aqueous 
solution. 

e The methanol that crosses over from the anode to the cathode 

is completely oxidized at the cathode catalyst layer. 

The effect of buoyancy in methanol transport is neglected. 

e The membrane (assumed to be Nafion® 117) is fully hydrated 
and is impermeable to gases. 

e The pressure gradient across the different cell layers is 
neglected. 

e Ohmic losses in gas diffusion layers, channels and bipolar 

plates are neglected. 


Some of this assumptions could be easily revised to incorpo- 
rate additional effects in future extensions of the model, but will 
be maintained here for simplicity. 


2.2. Physical domain 


We shall assume a parallel channel geometry for the anode 
current collector. Accordingly, when describing the flow in a 
single channel we shall use periodic boundary conditions at 
the channel/rib mid-planes to reduce the computational cost. 
The investigation of the effect of geometry — including cross- 
section channel shape and/or flow channel geometry — on fuel 
cell performance should be considered in future work. 

Fig. 1 shows an sketch of the physical domain considered 
here, which can be divided into seven regions: 


(1) anode channel (ac); 

(2) anode gas diffusion layer (agdl); 
(3) anode catalyst layer (acl); 

(4) polymer membrane (mem); 

(5) cathode catalyst layer (ccl); 

(6) cathode gas diffusion layer (cgdl); 
(7) cathode channel (cc). 


Since we are solving neither the electric field nor the tem- 
perature field, we omit from the computations both the anode 
current colector (acc) and the cathode current colector (ccc). The 
figure also shows the nomenclature used for the dimensions of 
the anode flow distributor (channel depth, ôac, channel width, 
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Fig. 1. Schematic of the modeling domains covered by the 1D and 3D models 
showing the coordinate system, and the notation used for the anode channel 
length, L, width, wac, depth, Sac, rib width wr, and the thickness of the different 
layers of the MEA. Top: side view; bottom: cross-sectional view. The 1D model 
does not account for the detailed geometry of the cathode flow distributor. 


Wac, rib width, wr) and the thickness of the different layers of 
the MEA. 


2.3. 3D model (anode channel and gas diffusion layer) 


The complete Navier-Stokes equations, complemented with 
the conservation equation for methanol, are solved in a 
three-dimensional (3D) domain with a control volume-based 
discretization to obtain the velocity and pressure distributions, 
as well as the concentration of methanol in the anode channel 
and gas diffusion layer. 


2.3.1. Flow field 

The Navier-Stokes and continuity equations governing the 
steady motion of an incompressible fluid of constant density 
and viscosity, p and u, through an isotropic porous medium 
may be written as [27]: 


V-u=0 (4) 


5u .Vju = —V p + Evu + Su (5) 


in terms of the porosity € of the porous matrix. Here we shall use 
the value éaga] = 0.6 for the porosity of the anode gas diffusion 
layer, while we simply set € = 1 in the anode channel, thus 


reducing (4) and (5) to the usual Navier-Stokes equations. In 
writing the above equations we have used the superficial velocity, 
u, based on the volumetric flow rate, and we have neglected 
buoyancy effects. Since we will incorporate the effect of the 
electrochemical reactions as boundary conditions at the anode 
gas diffusion layer/catalyst interface, we also assume that there 
are no mass sources or sinks associated with the electrochemical 
reactions inside the computational domain. 

As discussed in Section 2.1, we consider that the concentra- 
tion of methanol in the anode channel is so small that it does not 
alter significantly the physical properties of water. Accordingly, 
p and u represent the effective density and viscosity of water, 
given in Refs. [28,29] as 


p = 1000 — 0.0178(T — 277.15)!7 (kgm73) (6) 


u = 0.458509 — 5.30474 x 10777 + 2.31231 x 1075 T? 
—4.49161 x 1078 77+3.27681x107!! T+ (kgm7! s7!) 


(7) 


in terms of the temperature T of operation, which here denotes 
the absolute (Kelvin) temperature. 

The laminar flow in the porous gas diffusion layer is mod- 
eled by the addition of a momentum source in the governing 
momentum equation: 


sva-(H)s ® 


where K is the isotropic permeability of the medium. Note that 
for the typically small values of the permeability K of the porous 
regions of fuel cell electrodes (we use K = 107 !? m°) the addi- 
tion of this source term simply reduces (5) to Darcy’s Law. The 
source term, representing the extra resistance force suffered by 
the fluid due to the presence of the porous medium, vanishes at 
the anode channel. 


2.3.2. Methanol distribution 

Due to the small methanol concentration in the anode of a 
DMFC, the velocity field, u, is decoupled from the methanol 
concentration field, Cm, which is governed by the convec- 
tion/diffusion equation: 


V - (UCm) = DË V? Cm (9) 


where DSF is the effective diffusion coefficient of methanol in 
water, given by Yaws [30]: 


pet = 10754163—999.787/T m2 s7! 
pet z f m,ac (10) 


m pet 1.5 pet 


m,agdl — €agdl m,ac 
at the anode channel and gas diffusion layer, respectively. As can 
be seen, in the second expression we use the Bruggeman correc- 
tion factor to take into account the extra resistance to diffusion 
of the porous medium. 


2.3.3. Boundary conditions 
Eqs. (4), (5), and (9) must be integrated supplemented with 
appropriate conditions on the boundaries of the 3D domain 
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shown in Fig. 1. At the non-permeable channel walls we pre- 
scribe non-slip and zero gradient boundary conditions for the 
velocity and concentration fields, respectively: 


u=0, n-VCm=0 (11) 


where n is the outward normal vector and V = 0/dxe, + 
0/0ze, + 0/dze, is the gradient operator. At the inlet of the 
anode flow channel we specify the velocity and methanol con- 
centration: 


u = Uinex, Cm = Cmin (12) 


while at the symmetry planes and channel exit we set zero gradi- 
ent boundary conditions for both the velocity and concentration 
fields: 


n -Vu =Q, n-VCm=0 (13) 


At the porous portion of the inlet and outlet boundaries a zero 
gradient boundary condition is prescribed for both pressure and 
methanol concentration. The internal boundary condition at the 
anode channel/gas diffusion layer interface is appropriately han- 
dled by imposing the continuity of flux and scalar variables. 

The only non-trivial boundary condition is that imposed at the 
anode gas diffusion layer/catalyst interface, where we prescribe 
the molar flux of methanol: 


(u Cm — DRV Cm)| n= Nm (14) 
y= 


which is due to methanol crossover and to the methanol con- 
sumption by the anodic reaction, and the normal velocity: 


1\ iW 
aati (ns + 5) 1P (15) 


which is due to the electro-osmotic flux of water across the 
membrane and to the water consumption by the anodic reaction 
— both proportional to the local current density i generated at 
the anodic reaction. Here F is Faraday’s constant and Ww is the 
molecular weight of water. The electro-osmotic drag coefficient 
of water ny appearing in (15) is defined as the number of water 
molecules dragged by a hydrogen ion moving in the membrane, 
and is given by Guo and Ma [13]: 


n\ = 2.9 exp 1029 (= — =) (16) 

oa 333 T 
It is important to note that the local molar flux of methanol 
at the catalyst layer, Nm, and the local current density, i, that 
appear on the right hand side of Eqs. (14) and (15), are not 
known a priori. The purpose of the 1D model presented below 
is precisely to calculate the values of Nm and i at the active 
boundary from the analysis of the remaining elements of the 
cell, namely the membrane, the catalyst layers, and the cathode 
flow field, thereby closing the mathematical problem. 


2.4. 1D model (MEA and cathode) 


The 3D mathematical model for the anode flow field 
described above will be complemented with a one-dimensional 


(1D) across-the-membrane model that provides the local val- 
ues of the molar flux of methanol, Nm, and current density, i, 
at the active boundary in terms of the cell voltage, V, and the 
local methanol concentration there, Cm,aci = Cmly=0. As previ- 
ously discussed, this 1D model, including both catalyst layers, 
the membrane, and the cathode flow field, will allow us to close 
the mathematical problem through Eqs. (14) and (15). 


2.4.1. Transport of O2 in the cathode 

Unlike ordinary polymer exchange membrane fuel cells, 
DMEFCs are known to suffer from mass transport limitations 
mostly at the anode [3]. As a consequence, the anode flow field 
design turns out to be a key element for the optimal operation of 
DMEFCs [31]. In fact, with proper control of water crossover to 
avoid cathode flooding, oxygen transport in the cathode is typ- 
ically sufficient to sustain the — usually low — current densities 
generated in DMFCs, even in the case of air-breathing DMFC 
stacks [32]. This is especially true in the low current density 
regime, for which our single-phase model is expected to give 
the best results. 

To reduce the computational cost, we shall circumvent the 
analysis of the cathode flow field by introducing an overall mass- 
transfer coefficient to model the transport of O2 from the ambient 
air to the cathode catalyst layer. Thus, we simply write the molar 
flux of oxygen that reaches the cathode catalyst layer as 


No, = &2(Co3,amb — Co,ccl) (17) 
in terms of a mass-transfer coefficient 
1 8 ae 
ced 
a2 = | — + = (18) 
(= | 


that includes both the convective and diffusive transport losses 
through the mass transfer resistances ho, and DË cel /Scgdl, 
respectively. 

In the above expression, dcgqi is the thickness of the cathode 
gas diffusion layer and, using the Bruggeman correction factor, 
we write the effective diffusivity of oxygen in the cathode gas 
diffusion layer as 


DË cel = el D0» air (19) 
where €cgqi is the porosity of the porous matrix, and 
3/2 
T P, 
f amb 
Donas = Dla ag) ("8") (20) 


is the diffusivity of oxygen in air, expressed in terms of its 
reference value De air at 298 K and 1 atm, given in Table 1. 

It is worth noting that when imposing a constant value 
Co;,amb for the oxygen “ambient” concentration along the cath- 
ode channel/flow distributor, a large stoichiometric excess of air 
(or oxygen) is being implicitly assumed. To be consistent with 
this approximation, we shall also neglect the convective loses 
(ho, = oo) when evaluating the molar flux of oxygen reaching 
the cathode catalyst layer. Nevertheless, in real devices water 
flooding may obstruct oxygen transport in the cathode resulting 
in smaller values of ho, and DÖ egal (and therefore of œ2) than 
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Table 1 
Physical constants and design parameters involved in the 3D model for the anode 
channel (ac) + gas diffusion layer (agdl). 


Parameter Physical description Value Reference 
F Faraday’s constant 96485 C mol7! - 

Wm Molecular weight of methanol 0.032 kg mol! - 

Ww Molecular weight of water 0.018 kg mol! - 

bac Porosity of ac 1 - 

Eagdl Porosity of agdl 0.6 Assumed 
Kac Permeability of ac oo - 

Kagal Permeability of agdl 107!2 Assumed 
L Anode channel length 8 x 10-2 m Assumed 
Wac Anode channel width 1.5 x 10-7 m Assumed 
Sac Anode channel depth 1x 1073m Assumed 
Wy Anode land width 1x 107-3m Assumed 
Sagdl Thickness of the agdl 1.5 x 1074m Assumed 


those considered here. This could lead, mostly for high current 
densities, to mass transfer limitations at the cathode that are not 
included in this simplified model. 


2.4.2. Mass balance of O2 at the cathode catalyst layer 
The oxygen that reaches the cathode catalyst layer either com- 

bines with the electrons and protons, or reacts directly with the 
methanol that crosses the membrane, to form water according to 
reactions (2) and (3), respectively. Accordingly, the mass balance 
of O2 at the cathode catalyst layer can be written as 
N i 3 

O: = 4F + z Netoss (21) 


with Ncross and i given by Eqs. (23) and (27) below, respectively. 


2.4.3. Mass balance of methanol at the anode catalyst layer 

Similarly, the methanol that reaches the anode catalyst layer 
by convection and diffusion from the anode backing layer must 
equal the amount of methanol consumed at the catalyst layer 
by reaction (1) plus the molar flux of methanol that permeates 
through the membrane, according to 


i 
Nm = 6F ao Neross (22) 


with Necross and i given by Eqs. (23) and (30), respectively. 


2.4.4. Methanol crossover 
Methanol transport across the membrane is driven by molec- 
ular diffusion, pressure gradient and electro-osmotic drag [2]. 
However, the effect of pressure gradient is typically small and 
can be neglected in the first approximation [15]. Then, assuming 
Fickian diffusion for methanol in the membrane, the molar flux 
of methanol can be written as 
Neross = no — p% 2 


m,mem dy 


(23) 
mem 
where nī is the electroosmotic drag coefficient of methanol, 
defined as the number of methanol molecules dragged by a 
hydrogen ion moving in the membrane, and D ss is the effec- 


tive diffusion coefficient of methanol in the membrane (Nafion® 
117) assumed to be independent of the concentration, which is 


a good approximation for the low methanol concentrations that 
we shall consider here. 

The electroosmotic drag coefficient of methanol is given by 
Ren et al. [33]: 


W, 
ng = — ny Cm,acl (24) 


Pw 
in terms of the electro-osmotic drag coefficient of water, ny, 
defined in Eq. (16) above. 

Although there are different expressions available in the lit- 
erature for the effective diffusion coefficient of methanol in 
the membrane (see, e.g., the discussion given in Ref. [34]), we 
choose that given by Scott et al. [1]: 


1 1 
D hei = 4.9 x 107! exp [2436 (5 — z) (m? s7!) 


which roughly corresponds to the temperature range covered in 
this study, 298 K < T < 363K. 

Due to the fast oxidation of methanol in the cathode, the 
resulting methanol concentration at the cathode catalyst layer, 
Cm,ccl, is very small, and can be considered to be zero in the first 
approximation. Accordingly, the diffusive flux of methanol in 
(23) can be approximated as 


dCm ~ neff Cm,acl 


m,mem 5 
mem 


(26) 
mem 
where ômem is the thickness of the membrane. In a conven- 
tional cell dem is typically small compared to the channel width 
and length. Consequently, the methanol concentration gradients 
across the membrane (i.e. in the y-direction as defined in Fig. 1) 
are anticipated to be large compared to those in the x- and z- 
directions. Thus, the flux of methanol through the membrane is 
mainly determined by the local value of the methanol concen- 
tration at the anode catalyst layer Cm,ac1, as implicitly stated by 
Eqs. (23) and (26). 


2.4.5, Electrochemical reduction of O2 in the cathode 

We shall assume that the electrochemical reduction of oxygen 
at the cathode follows Tafel kinetics with first order dependence 
in oxygen concentration: 


ar | Coy cel ao F 

+i, =6 D ex 27 
Ut lp cel (Aig )e Coret p RT Nc (27) 
where 
ip = 6F Neross (28) 


is the parasitic current density due to the permeability of the 
membrane to methanol, i.e. the current that would be gener- 
ated by the methanol that crosses-over the membrane, ôccı is 
the thickness of the cathode catalyst layer, ac is the effective 
cathodic catalyst surface area per unit volume, Co, cci is the 
molar concentration of oxygen at the catalyst layer, œe is the 
cathode transfer coefficient, and ioc is the exchange current den- 
sity of the cathodic reaction, given as a function of T by Wang 
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and Wang [15]: 


. : 73200 / 1 1 
ine = it exp | 7 (5 =)| (29) 


in terms of its reference value iret at 353 K. 


2.4.6. Electrochemical oxidation of methanol in the anode 
In this study, we adopt the Meyers and Newman [19] kinetic 
model for the oxidation of methanol on Pt—Ru catalysts: 
K Cyacl CXP(@a F/ RTNa) 
Cmact + A exp(œa F/RTnha) 


i = daci(aio)a (30) 
where dac) is the thickness of the anode catalyst layer, aa is the 
effective anodic catalyst surface area per unit volume, iọ,a is the 
exchange current density of the anodic reaction, Cm,acı is the 
molar concentration of methanol at the catalyst layer, œa is the 
anode transfer coefficient, and « and à are two experimentally 
fitted coefficients [8]. 

In the above expression, the anode exchange current density 
is given by Wang and Wang [15]: 


l . 35570 / 1 1 
ina = i exp | 7 (= =) | GD 


in terms of its reference value e at 353 K, and the molar 
methanol concentration, Cy ac}, is assumed to be uniform across 
the catalyst layer, a condition that is generally satisfied due to 
the small thickness of the catalyst layer and the high rate of 
methanol transport across it [26]. 

The kinetic model (30) was derived by Meyers and Newman 
[19] assuming complete methanol oxidation. Similar expres- 
sions have been recently used by Garcia et al. [8] and Kulikovsky 
[26] as improved models for methanol electro-oxidation, a pro- 
cess that is known to deviate significatively from Tafel kinetics. 
In both cases, Eq. (30) was used in the context of fully 1D 
models. 

As previously discussed, Eq. (30) describes the transition 
from first-order kinetics for low methanol concentrations and 
high overpotentials to zeroth-order kinetics for high methanol 
concentrations and low overpotentials. According to Meyers and 
Newman [19], the reason for this transition is that the rate of elec- 
trochemical oxidation of methanol is mainly determined by the 
desorption of the reactive molecules from the catalyst surface, 
except for sufficiently low methanol concentrations (i.e. high 
current densities), when the diffusive transport to the catalyst 
surface becomes the rate-determining step. 

As a consequence, the above kinetic model can also be 
derived from a simple two-step reaction mechanism that 
accounts for a slow potential-independent step of methanol 
adsorption on the catalyst layer, coupled to a second step with 
Tafel kinetics corresponding to the electrochemical conversion 
of the adsorbed species [26]. 

It should be emphasized that the above kinetic expression 
(30) avoids the use of nonintuitive transitions between different 
reaction orders at certain threshold concentrations. In particu- 
lar, it shows that the threshold concentration depends on the 


anodic overpotential according to coe ~ à exp(a,Fn,/RT), 


and therefore cannot be taken as constant, as done in previous 
work [15]. 


2.4.7. Equation for the cell voltage 
The cell voltage V is determined by the equation: 
. mem 


V = Ecell — Na — Ne — i (32) 


Omem 


where Ee; is the ideal electromotive force of the cell, na and 
Nc are the anodic and cathodic overpotentials, respectively, and 
the last term represents the ohmic drop across the membrane, 
expressed here in terms of the local current density, i, the thick- 
ness of the membrane, dmem, and the ionic conductivity of the 
membrane, Omem, assumed to be a constant, since the membrane 
is fully hydrated in liquid-feed DMFCs. 

The ideal electromotive force of the cell is given by Scott et 


al. [1]: 
AN RT ( P ) 
log 
n F Pamb 
(33 
where E° 


cell İS the ideal electromotive force under standard con- 
ditions, i.e. Cm = 1000 mol m3, Pamb = l atm, T = 298 K, 
(dE/ IT )iig is the rate of change of Eeey with T, and the last 
term represents the effect of pressure in the cathode potential. 
We shall assume here that the methanol and water of the overall 
reaction are liquid, hence AN/n = —0.5/2. 

Finally, the membrane (Nafion® 117) conductivity is given 
by Scott et al. [1]: 


1 1 
ee ee 1268 (= - 3) (34) 


is the 


Ece = Ecel F (T To) 
oT liq 


in terms of the temperature T of operation, where oem 


reference ionic conductivity of the membrane at 298 K. 
3. Solution procedure 
3.1. Solution of the 1D model 


Although the nonlinear character of the equations involved 
in the 1D model precludes an analytical solution of the prob- 
lem, a numerical solution can be easily obtained using iterative 
methods. 

Assume that we are given the local values of Cm,acı and V at 
the anode catalyst layer. These will be provided at every iter- 
ation by the numerical solution of the 3D anode model. Then, 
appropriate manipulations of Eqs. (17), (21)-(23), (27), (28), 
and (30) yield the following set of analytic expressions for the 
unknowns i, Nc, Na, Necross, ip,» Noy, and Co, cc) as a function of 
Nm and Cyc: 


Nm — Cm acl DET vein /mem 
. N. A C r = 6F : 2 
i(Nm m,acl) 1+ 6 nF (Cm,acl) 
RT 6FNmCo, ref 
hia E in| 2: | (36) 
cUVm, Cmvacl oF Sccl (Ai e CO>,ccl( Nm, Cm,acl) 


770 M. Vera / Journal of Power Sources 171 (2007) 763—777 


Na(Nm, Cm,acl) 
RT Cm,acl i(Nm, Cim,acl) 


Qa F ôaci(aio)a K Cm,acl — À i(Nm, Cm,acl) 
i(Nm, Cm,acl) 
Neross(Nim; Cm,acl) = Nm = ae (38) 
6F 
ip(Nm, Cin,acl) = 6FNcross(Nm, Cm,acl) (39) 
3 
No,(Nm, Cm,acl) = 7 Nm (40) 
3 N, 
Coz,cct(Nm; Cmact) = Cozamb = 5 (41) 
a2 


where the electro-osmotic drag coefficient of methanol appear- 
ing in Eq. (35), ng (Cmacl)» is given as a function of Cm,aci by 
Eq. (24). 

Substituting Eqs. (35)-(37) in (32) provides the following 
non-linear relationship between Nm, Cm acl, and V: 


f(Nm, Cm,acl) = Eeen — V — na(Nm, Cm,acl) — Ne(Nm, Cm,acl) 


E T (42) 

Omem 
which can be readily solved for Nm using a Newton—Raphson 
method for given values of Cm,acı and V. After solving for Nm, 
the remaining unknowns can be obtained from Eqs. (35)-(41). 

For low methanol concentrations the denominator inside the 
logarithm of Eq. (37) approaches zero, yielding high anode over- 
potentials and complicating the numerical solution of Eq. (42). 
In this limit, we used as initial guess for Nm the asymptotic 
solution of Eq. (42) given in the Appendix, which proved to be 
a good technique to accelerate the convergence of the numerical 
method. 

Illustrative results from the numerical solution of the 1D 
model are presented in Fig. 2. The local current density, i, is 
plotted together with the anodic and cathodic overpotentials, na 
and nc, versus the local concentration of methanol at the active 
layer, Cm,acı, for different values of the cell voltage, V. As can 
be seen, for small methanol concentrations i grows linearly with 
Cmacl, until it saturates for Cea ~ À exp(aaFn,/RT) due to 
non-Tafel effects. For higher methanol concentrations the cur- 
rent density decreases first due to methanol crossover (more 
prominent for low current densities and high cell voltages) and 
later to mass transport limitations at the cathode — as reflected by 
the sudden rise of the cathodic overpotential — until it eventually 
vanishes. 


3.2. Numerical solution of the coupled 3D/1D model 


For a fixed value of the cell voltage V, the coupled 3D/1D 
model was solved using the commercial finite-volume-based 
CFD code FLUENT®6.2. The computational domain was 
discretized using a Cartesian-structured grid generated with 
GAMBIT® 2.0. The total number of volume elements was 
24,000 in the anode channel and 20,000 in the gas diffusion 


i (mA/cm?) 


C 


(mol/m3) 


m,acl 


Fig. 2. Variation of the local current density i (upper plot), anode overpotential 
Na, and cathode overpotential ne (lower plot) with the local concentration of 
methanol at the anode catalyst layer, Cm,aci, for different cell voltages V (solid 
lines). Vertical dashed lines corresponding to 0.1 and 1 M methanol concentra- 
tion and undiluted methanol are indicated for illustrative purposes. The dashed 
lines in the upper plot represent the asymptotic solution for small values of Cm,acl 
given in the Appendix. 


layer. To ensure a good resolution at the most critical regions, 
the grid points were clustered at the anode channel/gas diffusion 
layer interface and in the vicinity of the anode catalyst layer. 

The boundary conditions at the gas diffusion layer/catalyst 
interface given by Eqs. (14) and (15) were implemented in 
FLUENT® through the use of user defined functions (UDFs). 
This allowed us to solve the 1D model at every iteration using 
the local value of Cm at the active boundary to obtain the corre- 
sponding molar flux of methanol Nm. 

Once the numerical solution was converged, the average cur- 
rent density of the cell was calculated according to 


1 
I= ido (43) 
Aacl i: 


where the local current density i is given by Eq. (35) and the 
surface integral is extended over the whole surface area Aacı of 
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the anode catalyst layer. The average parasitic current density 
due to crossover, Ip, was calculated similarly. Then, another 
value of V was given and the process was repeated until the 
whole polarization curve was obtained. 

To check the accuracy of the numerical solution, a grid- 
refinement study was performed. The refined grid was obtained 
by halving the node spacing in each coordinate direction, 
increasing by a factor of eight the total number of grid elements. 
The numerical results obtained with the refined grid showed rel- 
ative errors smaller than 107° in the numerical evaluation of the 
average current density as compared to the coarser grid. Con- 
sidering the inherent limitations of the model (in particular, the 
uncertainties in the physical parameters) this level of accuracy 
was considered to be appropriate. Using the coarse grid, the con- 
vergence of a single point in the polarization curve took about 
30 min of CPU time using a Pentium IV processor at 2.67 GHz 
with 1 GB RAM. 


3.3. Physicochemical parameters 


In the development of the present model several parameters 
have been introduced, such as the effective diffusion coeffi- 
cients, gas diffusion layer porosities, protonic conductivity in the 
membrane, anode/cathode transfer coefficients, exchange cur- 
rent densities, etc. Most of them were taken from the literature, 
as summarized in Tables 1 and 2. 

As a final remark, it should be noted that some parameters 
do not enter individually in the model. Instead, they appear 
grouped together in products such as ôacı(aio)ak OF ôcci(aio)e- 
Accordingly, if we, for instance, increase the specific area while 


decreasing the thickness of the catalyst layers by the same 
amount, the results will remain unaltered. 


4. Results and discussion 
4.1. Case studies 


The different cases under study are summarized in Table 3. 
Simulations were performed for different methanol feed con- 
centrations, Cm,in, temperatures of operation, T, methanol flow 
rates, Uin, and cathode fluid—either ambient pressure air or pure 
oxygen. The reference case (Case 0) corresponds to the flow of 
a 0.5 M (500 mol/m?) methanol dilution at 1 mm/s inlet velocity 
and 80°C, using ambient pressure air as the cathode fluid. 


4.2. Overall cell performance 


Figs. 3-5 summarize the computed overall cell perfor- 
mance for the different cases under study, including polarization 
and power-—current curves, average parasitic current density 
due to crossover, and fuel utilization FU = 100 x I/(U + Ip)% 
— defined as the ratio of the electrochemically reacted fuel 
and total fuel used — as a function of cell current den- 
sity. Different figures illustrate the effect of methanol feed 
concentration (Fig. 3, cases 0-3, 12), the effect of tempera- 
ture (Fig. 4, cases 0, 4-7), and the effect of methanol flow 
rate (Fig. 5, cases 0, 8-11). Fig. 3 also shows the effect 
of using ambient pressure air or pure oxygen as cathode 
fluid. 


Table 2 
Physical constants, transport, kinetic and design parameters involved in the 1D model for the anode catalyst layer (acl) + membrane (mem) + cathode (ccl, cgdl, cc) 
Parameter Physical description Value Reference 
F Faraday’s constant 96485 C mol~! - 
R Universal gas constant 8.31Jmol~! K7! - 
Win Molecular weigth of methanol 0.032 kg mol~! - 
Ww Molecular weigth of water 0.018 kg mol~! - 
DE ait Diffusivity of O in air at 298 K 2.5 x 1075 m? s7! Perry et al. [35] 
X Oan O2 molar fraction in the cathode channel 0.21 (air) / 1 (pure O2) - 
COs. amb O2 molar concentration in the cathode channel X0,,amb Pamb/ RT - 
CO ret Reference O2 molar concentration 0.21 Po/RT mol m7 - 
Po Reference pressure for Co, „er 10° Pa - 
Pirih Ambient pressure 1.013 x 10° Pa - 
Oa Anode transfer coefficient 0.5 Murgia et al. [7] 
Oc Cathode transfer coefficient 1.2 Murgia et al. [7] 
da Anode catalyst surface area per unit volume 6 x104 m7! - 
dc Cathode catalyst surface area per unit volume 6 x104 m7! - 
K Experimental constant 7.5 x1074 Garcia et al. [8] 
À Experimental constant 2.8 x 1073 mol m~? Garcia et al. [8] 
int Anode exchange current density at 353 K 94.25 A m~? Ren et al. [33] 
igt Cathode exchange current density at 353 K 0.04222 A m~? Wang and Wang [15] 
Coen Membrane (Nafion) conductivity at 298 K 7.3Sm7! Scott et al. [1] 
2i Open circuit voltage at 298 K 1.213 V Wang and Wang [15] 
(GE/8T )iig Rate of change of E? with T -14x 10-4VK"! Scott et al. [1] 
Back Thickness of the acl 3x 10> m Assumed 
mem Thickness of the membrane 1.78 x 1074m Assumed 
Seel Thickness of the ccl 3x 107m Assumed 
Segdl Thickness of the cgdl 1.5 x 1074m Assumed 
Ecedl Porosity of cgdl 0.6 Assumed 
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Table 3 
Operational parameters for the different cases under study 
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Case no. Cmin (mol m~?) T(K) Uin (mm s~!) Cathode fluid 
0 (reference) 500 353 1 Air 
1 1000 353 1 Air 
2 250 353 1 Air 
3 125 353 1 Air 
4 500 363 1 Air 
5 500 343 1 Air 
6 500 333 1 Air 
7 500 298 1 Air 
8 500 353 0.5 Air 
9 500 353 2 Air 
10 500 353 4 Air 
11 500 353 8 Air 
12 500 353 1 <CE:BOLD>0</CE:BOLD>2 
As can be observed, the model reproduces appropriately (27) as 
the parametric trends reported in the literature, predicting an 
j i . RT X O2,pure RT 
increase in cell voltage (and therefore power output) for increas- Ane X log = log = 0.0395V 
Qc F XOo,air ack 0.21 


ing inlet methanol concentration, temperature, and methanol 
flow rate. As illustrated in the inset of Fig. 3, this trend is reversed 
at low current densities, when the cell voltage drop due to 
methanol crossover is larger for higher methanol concentration. 

Fig. 3 also shows that, within the limitations of the present 
model, the main effect of increasing oxygen concentration in 
the cathode by using pure oxygen (X0, pure = 1) instead of air 
(X0),air = 0.21) is to reduce the cathodic overpotential by a con- 
stant amount, which can be estimated approximately from Eq. 


0.8 


0.7 
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P 


I/ +I) (%) 


(44) 


at T = 353K, independently of the cell current density. This 
effect rises the overall cell voltage, therefore increasing the cell 
power output, but does not affect the limiting current density, 
which in our simplified model is mainly determined by mass 
transport limitations at the anode. 

The parasitic current due to cross-over predicted by the model 
shows also values and trends similar to those reported in the lit- 
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Fig. 3. Computed cell performance corresponding to T = 80°C, Uin = 1 mm s~!, and different methanol feed concentrations, Cmin = 0.125, 0.25, 0.5, and 1 M (as 
indicated): polarization curve (upper left), power—current curve (lower left), average parasitic current density due to crossover (upper right), and fuel utilization % 
(lower right) as a function of cell current density. Cathode fluid: ambient pressure air (—), pure oxygen (——) (shown only for Cm,in = 0.5 M). The solid dots in the 


lower right plot represent fuel utilization at maximum power. 
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Fig. 4. Computed cell performance corresponding to Cm,in = 0.5M, Uin = 1 mm s~!, and different temperatures of operation, T = 25, 60, 70, 80, and 90° C (as 
indicated). See caption of Fig. 3 for details. 
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Fig. 5. Computed cell performance corresponding to Cm,in = 0.5 M, T = 80°C, and different methanol feed velocities, Uin = 0.5, 1, 2,4, and 8 mm s7! (as indicated). 
See caption of Fig. 3 for details. 
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erature [8,15]. Thus, the leakage current density decreases by 
operating the fuel cell at low methanol concentrations and low 
temperatures, and decreases almost linearly with cell current 
density. The results also show a reduction in parasitic current for 
low methanol flow rates. It is worth noting that even at the mass 
transport limited current density the model predicts a notice- 
able (i.e. non-zero) cross-over current, which is always larger 
than approximately 5% of the limiting current density (so that 
FUmax < 95%). 

It is also interesting to note that the fuel utilization at maxi- 
mum power output, shown by solid dots in the lower right plot 
of each figure, reaches always values between 80 and 90%. 
These values, which increase with temperature and decrease 
with feeding methanol concentration, are also in agreement with 
previously reported results [25]. 


4.3. Methanol distribution in the anode 


3D models constitute an excellent tool for exploring the spa- 
tial distribution of reactants and current density in the fuel cell. 
Figs. 6 and 7 show the distribution of methanol for the refer- 
ence case at various cross-sections along the anode channel, 
at the channel symmetry plane, and at the anode gas diffu- 


0.5 


z (mm) z (mm) 


sion layer/catalyst interface, corresponding to three different cell 
voltages: V = 0.8, 0.4 and 0 V. As expected, due to mass trans- 
port limitations the region below the current collector shows low 
methanol concentrations. Moreover, the methanol concentration 
at the catalyst layer is seen to decrease noticeably when increas- 
ing the cell current density, approaching a small — but finite — 
value at the limiting cell current density. 


4.4. Anode flow field and water management 


All the flows simulated herein exhibit Reynolds numbers of 
order unity based on the channel depth, ô., and methanol feed 
velocity, U, e.g., for the reference case Re = pUôc/u ~ 2.7. 
The flow is therefore steady and hydrodynamically developed 
over the entire channel, except for a short entry region of char- 
acteristic length 6, « L. Moreover, the water lost through the 
gas diffusion layer/catalyst interface, due to both electroosmotic 
drag and water consumption at the anodic reaction, reduces the 
volumetric flow of water as it evolves downstream. 

For illustrative purposes, Fig. 8 shows the streamlines of the 
flow at the channel symmetry plane and porous backing layer. 
The streamlines, obtained numerically, correspond to the same 
operating conditions considered in the previous section, i.e. the 


0.5 
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Fig. 6. Methanol concentration contours (mol m~?) at various cross-sections along the anode channel (x = 10, 30, 50, and 70 mm, from left to right) corresponding 
to (a): V = 0.8 V, 1=0.4mAcm~, (b): V = 0.4 V, I = 85mAcm~?, and (c): V = OV, I = Imax = 119mAcm~?. Results obtained for fixed values of Cm,in = 
500 mol m~? (0.5 M), T = 80°C, and Uin = 1 mms~!. The shaded region represents the anode gas diffusion layer. 
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Fig. 7. Methanol concentration contours (mol m~?) at the channel symmetry 
plane (upper plot) and gas diffusion layer/catalyst interface (lower plot) corre- 
sponding to the cases V = 0.8, 0.4, and 0 V shown in Fig. 6(a)—(c), respectively. 
The shaded region represents the anode gas diffusion layer (upper plot), and the 
region below the rib (lower plot). 


reference case operating at V = 0.8, 0.4 and OV. In the plots, 
the shaded region represents the anode gas diffusion layer, while 
the thick dashed line shows the dividing streamline separating 
the fluid that exits the domain through the channel exit from that 
lost through the catalyst layer. 

As should be expected, the fraction of fluid that is lost through 
the catalyst layer grows for increasing values of the cell current 
density, Z, or decreasing cell voltages, V—in particular, the sim- 
ulations reveal that about 6% of the inflow is lost for V = 0.4 V, 
and up to 10% for V = OV, i.e. at the cell limiting current den- 
sity. For increasing values of the cell current density the diffusive 
transport of methanol through the porous layer becomes, there- 
fore, more and more assisted by the convective transport of water 


x (mm) 


Fig. 8. Streamlines of the flow at the channel symmetry plane corresponding 
to the cases V = 0.8, 0.4, and 0 V shown in Figs. 6(a)-(c), respectively. The 
shaded region represents the anode gas diffusion layer, and the thick dashed line 
shows the dividing streamline. 


towards the membrane, a mechanism that probably contributes 
to sustain higher limiting current densities. The numerical results 
reveal that the convective flow of water through the porous layer 
also contributes to the transverse transport of methanol to the 
region below the ribs, although this effect can not be seen in 
Fig. 8. 


4.5. Current density, methanol crossover, and fuel 
utilization distributions in the anode catalyst layer 


Fig. 9 shows the predicted distributions of current density, 
methanol crossover, and fuel utilization in the anode cata- 
lyst layer for the reference case, operating at V = 0.4 V, I = 
85mA/cm?. As a consequence of the uneven distribution of 
methanol at the catalyst layer, the reaction occurs mainly below 
the gas channel, while the region below the current collector 
remains relatively inactive due to the low values of the methanol 
concentration there. On the other hand, the high methanol con- 
centrations imposed near the inlet by the feeding stream induce 
high cross-over currents there, which locally lowers the fuel 
utilization. 


4.6. Experimental validation of the numerical model 


To validate the mathematical model we performed a series 
of simulations to compare the numerical results with the exper- 
imental data reported by Sundmacher et al. [18]. During the 
validation, the geometric parameters were varied from those 
in Table 1 so as to reproduce the experimental conditions in 
Ref. [18], namely: L = 3 x 107?m, wac = 2 x 1073m, w, = 
1 x 107° m, ôac = 2 x 107° m, ôagdi = Seay = 0.3 x 107° m, 
mem = 0.178 x 107? m. The operating conditions were also 
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Fig. 9. Local distributions of the current density, i, methanol flux due to 
crossover, ip, and fuel utilization, i/(i + ip), at the anode catalyst layer for 
the reference case, Cm,in = 500 mol m7? (0.5M), T = 80°C, Uin = 1mms7!, 
operating at V = 0.4 V, I = 85mAcm~?. 


chosen so as to reproduce those in Ref. [18], i.e. T = 343K, 
Uin = 0.54 mm/s, Pamb = 2.5 barg. The specific catalyst surface 
area, a = 2 x 10*m7!, and the kinetic parameter, à = 1.25 x 
1072 mol m3, were also modified from the values shown in 
Table 1 in order to fit the experimental data. 

Fig. 10 shows the experimental polarization curves given in 
Ref. [18] together with the numerical results provided by the 
mathematical model. As can be seen, the proposed model fits 
satisfactorily the experimental data in the low current density 
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Fig. 10. Comparison of the experimental results reported by Sundmacher et al. 
[18](symbols) and the results provided by the mathematical model (solid lines) 
under similar operating conditions, for different methanol feed concentrations: 
(CA) 0.125 M; (A) 0.25 M; (O) 0.5M; (0) 2M. 


regime, where the cell voltage is mainly determined by activation 
polarization losses. However, multi-phase flow at both the anode 
and the cathode, which is not accounted for in this simplified 
model, may be partially responsible for the lack of quantita- 
tive agreement between experimental and numerical results for 
the high current densities reached for large methanol concentra- 
tions (2 M). Note, however, that even for the highest methanol 
feed concentrations the numerical results are very similar to 
those reported by Xu et al. [20] in the same parallel channel 
configuration. 

As a final remark, it is worth noting that the model results 
slightly overpredict the cell voltage for near open circuit condi- 
tions. The reason is that we are considering irreversible kinetic 
laws for both the electro-oxidation of methanol at the anode, Eq. 
(30), and the reduction of oxygen at the cathode, Eq. (27), while, 
indeed, for sufficiently low current densities the reverse reactions 
start to play a role. Accordingly, we can not expect the model 
to predict accurately open circuit voltages. A detailed account 
on the modeling and operation of a DMFC under open-circuit 
conditions can be found elsewhere [36,37]. 


5. Conclusions 


A novel 3D/1D, isothermal, single-phase model has been 
developed for liquid-feed direct methanol fuel cells (DMFC). 
The mathematical model accounts for several physicochemi- 
cal processes that limit the performance of DMFCs, namely 
the 3D convective/diffusive transport of methanol in the anode 
channel and gas diffusion layer, the parasitic electro-oxidation 
of methanol at the cathode due to methanol crossover, and the 
complex non-Tafel kinetics of methanol electro-oxidation at the 
anode, simulated using the kinetic model introduced by Mey- 
ers and Newman [19]. Although this model has been previously 
dismissed due to its complex non-linear nature [7], we were 
able to implement it in a computationally efficient way by using 
appropriate initial guesses based on asymptotic expansions. 

When applied to the steady-state operation of a single cell 
DMFC with straight channels the model was able to predict 
methanol concentration distributions, polarization and power- 
current curves, local current density profiles, local and total 
water crossover, as well as effects of inlet temperature, methanol 
concentration and methanol flow rate. 

The mathematical model was solved using a commercial CFD 
package, leading to polarization curves for different methanol 
feed concentrations, temperatures and volumetric methanol flow 
rate that are in agreement with the experimental results found in 
the literature. 

Further model improvements, including multi-species trans- 
port, multi-phase flow, electronic and protonic potentials, and 
water and heat management are worth future investigation to 
improve the predictive capabilities of the complex physicochem- 
ical phenomena involved in liquid-feed DMFC. 
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Appendix. Asymptotic solution of the 1D model in the 
limit of low methanol concentrations 


The asymptotic solution to Eq. (42) in the limit Cm,acı > 0 
is given by 


aci (a i0)a K 
See. a 
jä i 6FA2 
. &a/Ap 
6FCo;,ref daci(A io)a K 4 D aah 
Secl(A 10)¢ Co ,amb 6FA dmem 
Qa F 
x Chai exp [SE a = v] (45) 
where 
f bacila ioa DY 
ymax _ | ac a m,mem | ç 
m 6FA + Smem m,acl 
ôaci(a io)a k Wwn 
+| “ Fh z d Chaa (46) 
W 


is the maximum molar flux of methanol that can be sustained 
for a given methanol concentration. The expression for N? 
was obtained by equating to zero the denominator of the frac- 
tion inside the logarithm in Eq. (37). Notice that to reach the 
maximum flux of methanol the cell may require the application 
of negative voltages, which is impossible in practical devices 
where V must be always positive. 

Once N*** is known, the deviation of the molar flux of 
methanol below its maximum value can be calculated by substi- 
tuting Nm = NMax — W in Eqs. (35)-(41) and expanding (42) 
as an asymptotic series for small values of W. The zeroth 
order term of the resulting expansion involves two logarithmic 
singularities, coming from (36) and (37), that must be coun- 
terbalanced by W. Imposing that the zeroth order term must 
vanish in order to kill the logarithmic singularities, we obtain 
Eq. (45). 

Finally from Eq. (35) we may write the local current density 
for small methanol concentrations as i ~ 6F Nm, with Nm given 
by Eq. (45). 
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